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ABSTRACT 

A  micron  n"*"  —  n  —  n"*"  silicon  diode  is  simulated  via  the  hydrodynamic  model  for  carrier 
transport.  The  numerical  algorithms  employed  are  for  the  non-steady  case,  and  a  limiting 
process  is  used  to  reach  steady  state.  The  novelty  of  our  simulation  hes  in  the  shock  capturing 
algorithms  employed,  and  indeed  shocks,  or  very  rapid  transition  regimes,  are  observed 
in  the  transient  case  for  the  coupled  system,  consisting  of  the  potential  equation  and  the 
conservation  equations  describing  charge,  momentum,  and  energy  transfer  for  the  electron 
carriers.  These  algorithms,  termed  essentially  non-oscillatory,  have  been  successfully  applied 
in  other  contests  to  models  the  flow  in  gas  dynamics,  magnetohydrodynamics  and  other 
physical  situations  involving  the  conservation  laws  of  fluid  mechanics.  The  method  here  is 
first  order  in  time,  but  the  use  of  small  time  steps  allows  for  good  accuracy.  Runge-Kutta 
methods  allow  one  to  achieve  higher  accuracy  in  time  if  desired.  The  spatial  accuracy  is  of 
high  order  in  regions  of  smoothness. 
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1.  INTRODUCTION 


A  micron  n'*'  —  n  — n"*"  silicon  diode  is  sirnnlated  via  the  hydro^Iynamic  model  for 
carrier  transport.  The  nnmerical  algorithms  employed  are  for  the  non-steady  case, 
and  a  limiting  process  is  used  to  reach  steady  state.  The  novelty  of  our  simulation 
lies  in  the  shock  capturing  algorithms  employed,  and  indeed  shocks,  or  very  rapid 
transition  regimes,  are  observed  in  the  transient  case  for  the  coiipled  system,  con¬ 
sisting  of  the  potential  equation  and  the  conservation  equations  describing  charge, 
momentum,  and  energy  transfer  for  the  electron  carriers.  These  algorithms,  termed 
essentially  non-oscillatoiy,  have  been  successfully  applied  in  other  contexts  to  model 
the  flow  in  gas  dynamics,  magnetohydrodynamics  and  other  physical  situations  in¬ 
volving  the  conservation  laws  of  fluid  mechanics.  The  method  here  is  first  order  in 
time,  but  the  use  of  small  time  steps  allows  for  good  accuracy.  Runge-Kutta  meth¬ 
ods  allow  one  to  achieve  higher  accuracy  in  time  if  desired.  The  spatial  accuracy  is 
of  high  order  in  regions  of  smoothness. 
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The  conservation  laws  for  charge,  momentum,  and  energy  for  each  carrier  can 
be  obtained  by  taking  the  first  three  moments  of  the  Boltzmann  equation  [see 
Blotekjaer  [3],  Rudan  and  Odeh  [23],  and  Cercignani  [4]].  This  will  result  in  a 
system  of  five  equations  with  fourteen  variables  in  three  dimensional  space.  These 
variables  are  defined  by  particle  density  n,  velocity  u,  internal  energy  e;,  heat  flux 
q,  and  symmetric  pressure  tensor  p.  In  order  to  close  the  system,  we  assume  that  p 
and  q  are  defined  via  (cf.  Gardner,  Jerome  and  Rose  [9]), 
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where  kf,  is  the  Boltzmann  constant,  m  is  the  mass  of  one  carrier,  and  T  is  the 
carrier  temperature  assumed  to  be  scalar.  The  moment  system  also  includes  aver¬ 
age  collision  terms,  denoted  C„,Cp,  and  Cu,,  respectively,  which  describe  the  rate 
of  change  of  mass,  momentum  and  energy.  They  will  account  for  recombination 
regeneration  processes,  electron-electron  and  electron-lattice  collisions,  and  transfer 
of  energy  between  electrons  and  oe.  Their  explicit  form,  for  the  application 
considered  here,  is  given  by 
Cn  =0, 

~P 

Cp  =  -  {P  :=  carrier  momentum  density) 

-{w-^nkbTo)  ^  ^ 

Cw  =  - - -  (w  :=  carrier  energy  density) 

Here  we  use  the  following  empirical  relations  [see  [9]]  for  relaxation  times  and  the 
thermoconductivity  constant: 


mpoTo 
eT  : 


^pokbTTo  ^ 

2evj{T  +  To)  2 


3po  kl  To 


K  = 


where  To  is  the  temperature  of  the  lattice  (possibly  a  function  of  space),  /io  is  the 
electron  mobility,  e  the  charge  of  an  electron,  and  Vg  is  the  saturation  velocity.  We 
shall  discuss  the  foundations  of  these  relations  below.  We  note  here,  for  consistency, 
that  e/  is  proportional  to  T,  and  in  the  following  section  we  shall  define  P  and  w, 
and  describe  the  equations  of  the  hydrodynamic  system  and  the  system  boundary 
conditions. 

The  recent  literature,  pertaining  to  the  device  model  simulated  here,  includes, 
in  addition  to  [23]  and  [9],  the  studies  by  Baccarani  and  Wardeman  [1],  by  Odeh, 
Rudan,  and  White  [19],  and  by  Fischetti  and  Rudan  [8].  In  each  case,  velocity 
overshoot  effects  are  simulated,  which  cannot  be  recovered  from  the  standard  drift- 
diffusion  model.  As  is  now  widely  understood,  the  effect  is  strongly  dependent  on 
the  semiconductor,  and  is  much  more  pronounced  in  gallium  arsenide  than  siUcon, 
for  example.  Since  velocity  overshoot  is  related  to  variations  in  termperature  and 
electric  field,  the  energy  equation  in  the  hydrodynamic  model  thus  allows  for  the 
heat  flux  term  q  governed  by  the  Wiedemann-Franz  law  [2],  (note  that  K.n  is  the 
conductivity  in  this  law)  and  the  coupling  to  the  potential  equation  is  essential  in 
this  case.  Moreover,  the  momentum  equation  plays  the  role  of  a  refined  constitutive 
current  relation.  Rather  than  seeking  modifications  of  the  classical  current  relation 
by  the  incorporation  of  relevant  physical  effects  beyond  critical  values  of  parameters 
such  as  k  ^  (see  [1]  for  such  values  and  Thomber  [30]  for  a  sophisticated  discus¬ 
sion  of  such  augmented  relations),  we  find  here  that  current  is  induced  by  variable 
concentration  and  momentum  within  the  system. 
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The  validity  of  the  assumptions  concerning  the  collision  terms  Cp  and  Ca,, 
termed  phenomenological  approach  in  [7],  has  been  discussed  by  Nougier  et  al.  [18], 
where  agreement  within  10%,  for  n  type  silicon  and  p  type  germanium,  has  been 
demonstrated,  relative  to  a  direct  simulation  of  the  Boltzmann  equation,  for  a  wide 
range  of  electric  fields  and  lattice  temperatures.  Further  assumptions  employed 
in  the  model,  also  used  in  [1],  [23],  [19],  and  [9],  are  the  existence  of  a  scalar 
quantity  for  the  temperature  tensor,  constructed  from  the  random  velocity,  and 
the  parabolic  representation  for  the  energy,  in  terms  of  its  kinetic  and  thermal 
components.  Following  Hess  and  Sah  [13],  and  [1],  we  have  assumed  a  constant 
diffusivity  coefficient,  which  is  consistent  with  the  assumptions  of  constant  effective 
mass,  and  scalar  temperature  and  diffusivity.  This  has  the  effect,  via  the  Einstein 
relation,  of  describing  the  electron  mobility  as  an  inverse  temperature  dependent 
function. 

Thus,  in  this  approach,  mobility  does  not  directly  depend  upon  velocity.  This 
mobility  relation  explicitly  dictates  the  expression  for  the  temperature  dependent 
momentum  relaxation  times  expressed  above,  as  well  as  the  choice,  exponent  = 
—  1,  in  the  Wiedemann-Franz  law  for  the  thermal  conductivity.  The  temperature 
dependent  energy  relaxation  time  is  deduced  by  using  the  empirical  expression  for 
the  mobility,  presented  by  Caughey  and  Thomas  [5],  so  aa  to  express  drift  velocity 
in  terms  of  temperature.  The  resulting  expression  is  therefore  empirical.  The 
steady-state  character  of  both  the  resulting  relaxation  time  expressions  has  been 
emphasized  in  [l].  Since  we  simultaneously  obtain  information  for  the  transient 
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analysis,  the  latter  results  must  be  interpreted  provisionally  with  respect  to  Tp  and 

Tu,- 

The  earliest  theoretical  predictions  of  velocity  overshoot  appear  to  be  due  to 
Rees  [20],  for  germanium,  and  to  Ruch  [22],  for  silicon  and  gallium  arsenide.  Fur¬ 
ther  studies  by  Shur  [26]  and  Maloney  and  Frey  [16]  gave  a  detailed  qualitative 
understanding  of  the  phenomenon;  in  particular,  it  was  understood  that  the  physi¬ 
cal  principle  underlyii.g  velocity  overshoot  is  the  finite  time  taken  before  the  energy 
of  the  carriers  reaches  equilibrium  with  the  new  field,  as  they  traverse  a  region  of 
rapidly  varying  electric  field.  These  ideas,  and  the  relation  to  ballistic  transport, 
have  been  discussed  and  critiqued  by  several  authors  (see  Shur  and  Eastman  [27] 
and  Hess  [12]). 

Increasing  sophistication  of  the  model  followed.  While  [26]  and  [16]  were  based 
upon  the  quasi-diode,  by  1982  Cook  and  Frey  [6]  had  incorporated  many  essential 
features  of  the  hydrodynamic  model,  coupled  to  the  potential  equation,  in  their 
simulations.  However,  they  did  make  the  assumptions  not  made  here,  of  vanishing 
heat  flux  and  negligible  convective  energy  term.  The  former  assumption,  reducing 
order,  substantially  alters  the  mathematical  character  of  the  energy  equation  in  ad¬ 
dition  to  its  physical  implications;  these  mathematical  implications  are  manifested 
principally  in  two  or  more  dimensions,  however.  Other  attempts  to  incorporate 
energy  and/or  transport  effects  followed,  based  upon  various  models.  Thus,  the 
foundation  for  the  work  of  Tang  [29]  was  the  spectral  method  introduced  by  Strat- 
ten  [28],  where  a  spherical  harmonic  expansion  of  the  distribution  function  leads 


5 


to  momentum  and  energy  relations  similar  to  those  subsequently  derived  for  the 
hydrodynamic  model..  Also,  McAndrew,  Singhal,  and  Heasell  [17]  employ  a  model 
involving  the  momentum  equation,  but  not  the  energy  equation. 

The  models  employed  in  [1],  [23],  19],  and  [9]  are  virtually  the  same,  with  some 
variations  due  to  mobility- temperature  representations,  and  the  joint  connection 
with  relaxation  times  and  the  Wiedemann-Franz  law.  The  model  in  this  paper  is 
the  same  as  that  of  [9].  The  numerical  methods  in  this  paper,  however,  are  quite 
different  from  those  in  previously  reported  work.  Both  the  parametric  regimes  and 
the  iterative  differencing  schemes  used  in  previous  work  were  chiefly  keyed  to  the  el¬ 
liptic  or  elliptic /parabolic  classification  of  the  underlying  equations,  in  other  words, 
to  what  might  loosely  be  characterized  as  the  subsonic  transport  regime,  where  con¬ 
centrations  and  velocities  would  not  be  expected  to  experience  shock  discontinuities. 
In  steady-state,  the  relevant  “soundspeed”  is  y/KsT/m,  obtained  by  linearization  of 
the  density-momentum  subsystem.  If  the  density- momentum- energy  subsystem  is 
linearized,  the  traditional  “soundspeed”  of  is  obtained,  but  this  involves 

the  parabolic  mode  associated  with  the  energy  equation.  The  analysis  of  [9],  in  fact, 
shows  that  a  damped  Newton  method  can  be  justified  by  energy  type  estimates  in 
the  subsonic  case.  The  shock  capturing  methods  used  herein  have  both  the  ad¬ 
vantage  of  simulating  broad  parameter  ranges,  including  transonic  and  supersonic 
flow,  as  well  as  being  readily  extendible  to  two  dimensions.  Since,  as  noted  by  [1], 
velocity  overshoot  is  maximally  exploited  when  cold  carriers  are  injected  at  high 
field,  the  simulation  should  have  the  capability  of  covering  a  broad  temperature 
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The  above  scheme  can  be  extended  from  scalar  equations  to  systems  of  con¬ 
servation  laws  in  many  ways.  We  used  the  following  version  in  our  calculations  [see 
[25]]. 

Let  Ut  +  f{u)x  =  0  be  the  system.  Then  let  represent  the  Jacobian  matrix 
of  the  system. 

The  above  system  has  a  complete  set  of  eigenvectors  and  there  is  a  matrix  U 
such  that: 

U~^^U  =  A 

ou 

where  A  is  a  diagonal  matrix. 

Let  =  ^3^  and  6  =  ^  and  B  =  ^  and  h  =  -I-  jv^ .  Here,  c  is  the 

velocity  of  sound  waves  in  the  electrons  and  in  this  connection  we  define  the  Mach 
number  of  the  fluid  as  M  =  ^.  For  our  system  the  above  matrices  are  explicitly: 


/I  1  1  \ 

U  =  I  m(u  —  c)  mv  m(v  -f  c)  1 
\  m(—cv  -i-  h)  m(ct>  4-  k)  J 


U-\=  l-B 
\  hiB  - 


-If  J_  j. 

2  V  cm  '  m  ' 
vb 
m 

i(  j_  _ 

2  '  cm  m  ' 


/v-c  0  0  \ 

A  =  (A, 3)  =  I  0  V  0  ) 

\  0  0  V  +  c  j 


In  order  to  do  the  field  by  field  decomposition  of  the  fluxes  we  evaluate  the  U 
matrix  at  the  Roe  average  of  and  Uj+i  [see  [21]]. 
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Let  then  we  calculate 

This  gives  us  the  components  of  the  divided  differences  of  the  f  in  the  direction 
of  the  fields; then  using  the  algorithm  for  scalar  equation  we  calculate  the  flux  for 
each  field.  We  use  =  A(u^^i)  for  each  field.  Then  we  transform  back  to  the 
original  variables  via: 

This  scheme  will  be  r-th  order  accurate  in  space  where  r  is  the  order  of  the 
reconstructed  polynomial  Rj^^{x).  To  obtain  high  order  accuracy  in  time,  having 
the  ENO  property,  one  may  use  the  Runge-Kutta  methods  developed  in  [25].  How¬ 
ever,  since  our  primary  interest  here  was  in  steady  state  calculations  we  used  only 
first  order  forward  differencing  in  time.  The  resulting  scheme  would  be  unstable  if 
the  stencils  were  fixed;  the  adaptivity  stabilizes  the  linearly  unstable  method. 

4.  THERMAL  EQUILIBRIUM 

If  we  fix  the  bias  voltage  to  be  zero  we  expect  the  system  to  reprh  an  equilibrium 
state  where  the  temperature  of  the  electron  converges  to  the  lattice  temperature 
and  velocity  to  zero. 

Under  these  assumptions  the  system  will  be  simplified  to 

khTrij:  =  en$j. 

=  e(n  -  no) 
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TlriO)  =  nx(l)  =  0 


e  rtf  e  n. 


Note  that  we  get  the  following  relation: 


r~^y  ^ ) 

—  n  ^  k  t,T  ' 


n{x)  —  n,e^ 


If  we  let  'I'(x)  =  and  7^  =  ^7^  then  the  above  system  is  equivalent  to 

the  following  non-linear  boundary  value  problems: 


I  -  Hic'*'-')  =  -no(x) 

I  '^,(0)  =  ^,(1)=0 

(  ^'^ln{n(x))ii  -  n{x)  =  -nol-O 

I  n,(0)  =  n,(l)=0 

Rather  than  solving  any  of  the  above  we  choose  the  following  initial  conditions  and 
let  the  system  evolve  in  time. 

n(x,0)  =  noi-r)  e(x,0)  =  0  r(x.O)  =  ro 


The  calculated  solution  contains  the  junction  and  its  electric  field.  Also,  the 
density  has  to  be  used  as  initial  condition  for  the  rest  of  the  calculations.  The 
numerical  results  are  presented  in  the  numerical  section. 

5.  PHYSICAL  PARAiMETERS  OF  THE  DEVICE 

The  simvilated  device  is  a  silicon  MOSFET  channel.  The  device  is  an  —  n  — 
junction  of  length  one  micron. 

We  use  the  following  doping  profile: 

)  =  5  X  10' '  i  f  0  <  x  <  .25  or  .75  <  x  <  1 
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n£)(x)  =  2  X  lO^^cm  ^  if  .35  <  x  <  .65 

The  two  regions  are  connected  by  a  properly  scaled  version  of  the  following  poly¬ 
nomial; 

Q{x)  =  — 5x‘  -|-  21x^  —  35a;^  -f  35x  +  16. 

We  also  use  the  following  conditions  at  the  boundary: 

kkT  .n(0) 

- ln{ - ) 

e  Tii 

$(1)  =  — — ln{  ^  ^ +  vbtas{t) 
e  n, 

nj.(0)  =  ni(l)  =  0  p^(0)  =  Pj:(1)  =  0  u’j.(0)  =  1 )  =  0 


We  use  the  following  system  of  units: 

Length  :  10“®m  =  Imicron  Time  :  =  Ipico  second 

Mass  :  Ab;  Potential  :  Volt  Temperature  :  Kelvin 

Energy  :  10~'**J  Charge  :  10~'®C  Capacitance  ; 

In  the  above  system  of  units  the  constants  have  the  following  values; 

m  —  0.26me  =  .26  x  .9109  po  =  0.14  e  =  .1602  r,,  =  .1 

A:;,  =  .138046  X  10~^  n,  =  .014  e  =  11.7  x  S.S5418 

6.  NUMERICAL  RESULTS 

In  the  hrst  experiment  we  calculate  the  density  of  electrons  under  zero  l)ias.  The 
result  can  be  used  for  I’alculating  the  junction  capacitance  and  voltage  diffi'n'iice 
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across  the  junction.  Also  we  use  the  calculated  density  as  the  initial  data  for  the 
rest  of  the  experiments. 

From  a  numerical  point  of  view  this  experiment  can  be  used  to  test  the  accuracy 
of  the  scheme.  Since  the  velocity  and  current  must  go  to  zero,  the  difference  can  be 
used  as  a  measure  of  error. 

In  our  experiment  velocity  and  current  converged  to  a  small  oscillatory  solution 
around  the  zero  value.  Also  we  expect  the  oscillations  in  the  velocity  field  to  create 
some  difficulty  for  our  numerical  scheme  since  it  could  lead  to  differencing  in  the 
wrong  direction. 

Indeed,  this  created  spurious  oscillations  that  were  eliminated  by  a  small  mod¬ 
ification  in  the  ENO  scheme.  We  restricted  the  movement  of  the  stencil  such  that  it 
would  interpolate  across  oscillations.  This  can  be  done  in  the  following  fashion.  For 
the  sixth  order  ENO  note  that  the  ishift  variable  takes  values  between  -5  and  1.  We 
simply  forced  it  to  be  between  -4  and  0.  This  procedure  would  create  small  oscilla¬ 
tions  across  shocks,  but  one  can  use  a  modified  version  to  overcome  that  difficulty. 
In  that  version  for  updating  ishift  we  use; 

‘  \rij  +  ishift)!  >  |((1  +  e)f'{j  -I-  ishift  —  1)|  then  ishift=ishift  -  1, 

vV'l-’ere  e  is  of  order  of  the  truncation  error  of  the  scheme  and  positive  or  negative 
depending  on  the  value  of  ishift.  This  wouH  make  the  choice  of  stencil  biased 
towards  the  center.  In  numerical  calculations  this  modification  led  to  a  smoother 
solution  with  fewer  spurious  oscillations. 
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The  results  from  the  the  experiment  are  shown  in  the  fig  2-6.  In  figure  1  the 
density  of  donors  is  shown  for  reference.  As  can  be  seen  in  figure  2  ,the  electrons 
from  both  sides  expand  into  the  channel.  This  creates  electron  waves  which  settle 
down  approximately  in  three  pico  seconds.  The  junction  is  formed  and  the  electric 
field  over  the  junction  is  created.  In  figure  4  the  steady  state  of  the  velocity  is 
shown.  The  amplitude  of  the  oscillation  indicates  the  size  of  the  error  in  the  rest  of 
the  calculations. 

In  the  next  set  of  experiments  we  apply  one  volt  at  time  zero  and  calculate  the 
evolution  of  the  system  in  time.  The  system  reaches  steady  state  in  approximately 
five  pico  seconds.  In  figure  7  the  current  is  shown  as  a  function  of  time,  and  in 
figures  8  and  9  the  evolution  of  electric  field  and  velocity  are  shown.  In  figures  10 
and  11,  for  comparison,  we  duplicat'^  the  results  obtained  by  [19],  [9]. 

As  an  approximation  to  an  TV  curve  we  apply  a  bias  voltage  in  form  of  a 
“ramp”.  Then  we  graph  the  calculated  current  versus  voltage  in  figure  12.  If  th< 
slope  of  the  ramp  is  small  it  would  be  very  close  to  an  TV  curve.  There  is  some 
oscillation  in  the  curve  which  is  due  to  the  jump  in  the  derivative  of  the  bias  voltage. 
If  the  slope  of  the  ramp  is  small,  it  would  disappear.  In  figures  13-17  we  show  the 
evolution  of  the  solution  as  the  velocity  is  increased. 

If  the  temperature  of  the  lattice  is  raised  the  sound  speed  increases  and  the 
Mach  number  becomes  small.  To  test  the  scheme  for  small  Mach  number  we  set 
the  lattice  temperature  at  lOOOK.  The  evolution  of  velocity  and  temperature  are 
shown  in  the  next  two  graphs. 
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In  the  next  experiment  we  set  the  lattice  temperature  to  50K.  This  would 
reduce  the  sound  speed  .  The  results  are  given  in  figures  20-22.  We  emphasize  that 
we  need  to  modify  our  choice  of  stencil  to  avoid  oscillations  across  steep  gradients 
or  shocks. 

In  the  figures  23  and  24  we  show  just  such  a  situation  where  the  solution 
develops  a  very  sharp  profile.  The  density  and  the  velocity  are  shown  here.  The 
shock  capturing  algorithm  is  designed  to  handle  such  situations. 

In  the  next  experiment  we  set  our  bias  voltage  to  be  a  square  wave  and  then  a 
sinusoidal  wave.  The  results  are  shown  in  figures  25  -29. 

In  the  last  experiment  we  compare  the  calculated  solutions  in  the  steady  state 
using  different  orders  of  accuracy  in  the  x  variable.  From  extensive  calculations  we 
found  out  that  momentum  is  the  most  sensitive  variable  to  error  and  comparison 
of  momentum  would  show  errors  that  are  not  obvious  in  the  other  variables. 

The  calculated  results  show  spurious  oscillations  for  low  order  schemes  that 
(perhaps  surprisingly)  decrease,  with  increasznr/ order  of  accuracy. 

7.  CONCLUSION 

We  presented  a  numerical  scheme  that  is  able  to  solve  the  Hydrodynamic  model 
efficiently  and  accurately. 

The  .scheme  was  first  order  in  time,  but  it  can  be  extended  to  high  orders  of 
accuracy  using  the  Runge-Kutta  methods  developed  in  [25].  Also  the  high  order 
accuracy  in  space  turned  out  to  be  essential  at  the  junctions.  The  equations  allow 


very  steep  gradients  including  discontinuous  solutions,  but  there  was  no  observation 
of  steady  state  shocks.  However,  if  a  large  bias  voltage  is  applied  suddenly,  shocks 
form  which  dissappear  in  a  short  time. 

There  is  current  work  imder  progress  to  extend  the  calculations  to  holes  and 
also  to  two  spatial  dimensions. 

FIGURES 

Fig  1:  density  of  donors 

Fig  2:  velocity 

zero  bias.  Time  =  3,  ^  =  .2,  Aar  =  .01 
Fig  3:  temperature 

zero  bias.  Time  =  3,  ^  =  .2,  Ax  =  .01 
Fig  4:  velocity  in  steady  state 

zero  bias.  Time  =  3,  ^  =  .2,  Ax  =  .01 
Fig  5;  n  —  n£)  in  steady  state 

zero  bias,  T  =  5,  ^  =  .2,  Ax  =  .01 
Fig  6:  electric  field  in  steady  state 

zero  bias,  Time  —  3,  ^  =  .2,  Ax  =  .01 
Fig  7:  Current  vs.  Time 

v(t)  =  H(t),  Time  =  5,  ^  =  .2,  Ax  =  .01 
Fig  8:  electric  field  in  time 

v{t}  =  H(t),  Time  =  5,  ^  =  .2,  Ax  =  .01 
Fig  9:  velocity  in  time 
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range,  which  suggests  regions  of  transonic  flow  and  complex  behavior. 

2.  SYSTEM  AND  BOUNDARY  CONDITIONS 

The  moment  equations  are  specified,  via  the  summation  convention  in  three 
dimensions,  as 

5t(n)  +  ai.(nu,)  =  C„ 

dt{noj)  -i-  dx.{nviVj  +  Ptj)  =  nFj  +  Cj  component 

+  e/))  +  +  ei)  +  p,jVj  +  q,)  =  nF,v,  +  C^, 

where  all  of  the  terms  have  been  defined  except  the  forcing  function  F  due  to  the 
electric  field. 

The  Poisson  equation  for  the  electrostatic  potential,  the  first  of  the  Maxwell 
equations,  is  given  in  a  one  carrier,  electron  system  by 


=  e(n  —  rif?) 


where  is  the  density  of  donors.  Thus  the  forcing  function  F  becomes: 

F  =  —E  =  — V't, 
m  m 

We  define  P  =  mnv  and  w  =  mnej  +  =  jnkbT  +  jmnv^.  Let  the  vector 

u  =  (n,  P,  to)  denote,  respectively,  density,  momentum,  and  energy  of  electrons.  We 
get  the  final  form  of  our  equations  for  electrons  and  electric  potential  as: 

+  dxiinvi)  =  Cn 
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dt{Pj)  +  dxi(PiVj  +  hnT)  =  ~enEj  +  Cp 


dt{w)  +  dc.{vi{w  +  kbuT))  =  —enviEi  +  Cu,  4-  di-iundi.T) 


eV^<5  =  e(n  —  no) 


We  can  write  the  equations  in  the  compact  form  in  {x,y,z)  space  as: 


y-t  +  fi{u)x  +  f2{u)y  +  h{u)x  =  C(ii)  +  G(u,$)  +  (0,0,  V(/cnVT)) 


=  e(n  —  no) 


In  this  paper  then  we  concentrate  on  solving  the  following  system; 


=e(n  -  no),  E  = 

Tit  +  {nv)i  =  0, 


Pt  +  {Pv  +  k(,nT)j.  =  -enE - , 


Wt  +  {vw  +  vkbnT)i  =  —envE  — 


w  -  ^nkbTo  ^  ^  ^ 

+  (  KTll^  ) j 


with  the  following  boundary  conditions 


$(0)  = 

e  n, 

$(1)  =  — — ln(  -  4-  vbias{t) 

e  rii 

Wr(0)  =  Ui(l)  =  0. 


This  is  an  incompletely  parabolic  system,  very  similar  to  the  Navier-Stokes 
equations  of  compressible  flow.  Boundary  conditions  of  this  type  for  that  system 
have  been  discussed,  e.g.  in  [10]. 
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3.  NUMERICAL  SCHEME 


The  numerical  scheme  is  based  on  the  essentially  non-oscillatory  shock  captur¬ 
ing  algorithms  developed  in  [11].  We  use  a  sixth  order  stencil  in  space  and  forward 

u"— u" 

Euler  in  time.  Let  u"  =  u(jAx^nAt)  =  u(xj,tn)-  Also  let  D-u'J  =  ->  , 

n  ,,n  _  “7+1-“"  o  ..n  _  “7+1 -“7-1 

~  Ax  -  2 Ax  • 


Then  using  the  above  conventions  we  can  write  our  scheme  as: 


(2a) 


(2b) 


(2c) 


=  e(n"  -  no, ) 


EJ  = 


=  u?  - 


At 


— (/,+  ^  -  /,_ . )  +  AtC]  -H  AtG]  +  AtD.{nD+T) 


We  calculate  the  potential  from  the  density  and  then  we  calculate  the  electric 
field.  We  calculate  the  solution  in  the  next  time  step  from  the  electric  field  and  the 
’’flux”  calculated  in  the  previous  step. 

The  “flux”  term,  /'*_  i  ,  represents  the  numerical  flux  flowing  into  the  jth 
computational  cell  [see  Lax  [15]]. 

The  finite  difference  scheme  (2c)  is  an  approximation  of  ( lb)  .  The  system  ( lb) 
is  a  hyperbolic  system  of  conservation  laws  with  a  right  hand  side  which  consists 
of  a  viscosity  term  in  the  third  equation,  a  lower  order  term  and  a  non-local  term 
which  comes  from  solving  the  elliptic  equation  for  $  in  terms  of  n  and  ho  [see  9]. 

It  is  well  known  that  solutions  of  a  system  of  hyperbolic  conservation  laws  de¬ 
velop  discontinuities  after  a  finite  time  for  general  initial  data.  These  discontinuties 
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represent  shocks  and  contact  discontinuities.  In  the  present  problem,  these  may 
be  slightly  smeared  because  of  the  (extremely  small)  viscous  terms. There  are  ad¬ 
ditional  complications  due  to  the  stiffness  of  the  source  terms.  Also  these  viscous 
shocks  can  only  appear  if  the  Mach  number  exceeds  one. 

There  has  been  a  large  amount  of  work  in  recent  years  towards  developing  high 
resolution  approximations  of  the  numerical  solution  to  these  problems.  Scheme?  of 
the  general  type  (2)  are  called  shock  capturing  algorithms,  since  no  special  treatment 
of  the  discontinuties  is  needed. 

First  order  accurate  difference  schemes  are  generally  unsuitable  for  systems 
such  as  (lb)  due  to  the  extreme  smearing  of  the  steep  gradients  in  the  solution. 
Classical  second  order  accurate  schemes,  however,  are  generally  unsuitable  owing 
to  the  spurious  oscillations  which  they  generate  in  the  neighborhood  of  discontinu¬ 
ities.  Such  oscillations  are  not  only  aesthetically  undesirable  but  often  lead  to  a 
breakdown  in  the  stability  of  the  scheme  and  nonphysical  effects  such  as  negative 
densities. 

The  feature  of  the  classical  schemes  which  leads  to  unwanted  oscillations  is  that 
they  are  based  on  a  constant  stencil.  Moreover  these  schemes  are  linear  whenever 
the  equations  they  approximate  are  linear.  In  the  last  fifteen  years  there  has  been 
much  work  on  overcoming  such  problems  where  the  coefficients  and  the  width  of 
the  stencil  depend  on  the  current  solution  at  each  time  step.  Thus  the  schemes  are 
nonlinear  even  when  the  equations  are  linear.  This  is  absolutely  necessary  [11]  for 
the  properties  of;  both  high  order  accuracy  in  the  region  of  smoothness  and  the 


10 


absence  of  spurious  overshoot  and  undershoots  near  discontinuities,  to  be  valid. 

The  beisic  method  we  shall  use  here  which  has  these  desirable  properties  is 
a  modification  of  the  essentially  non-oscillatory  (ENO)  schemes  developed  in  [11]. 
The  key  ideas  are  the  following: 

In  each  computational  cell  the  data  at  t  =  <"  are  approximations  to  the  cell 
averages  of  the  vector  u,  denoted  by  u".  From  these  cell  averages  a  polynomial 
approximation  to  u,  of  order  r,  is  reconstructed  in  a  nonlinear,  adaptive  fashion. 
This  will  give  iis  an  r-th  order  accurate  method.  Values  of  for  s  near  zero 

are  used  in  such  a  way  as  to  minimize  the  likelihood  that  the  stencil  crosses  over  a 
discontinuity.  Finally  the  solution  to  this  piecewise  polynomial  initial  value  problem 
is  computed  approximately,  keeping  its  non-oscillatory  properties,  and  for  = 
t"  +  At,  At  small,  the  solution  is  reaveraged  over  the  cell  yielding  . 

As  an  example,  the  approximate  solution  to  the  scalar  equation  uj  =  — u,.. 
using  a  second  order  accurate  ENO  scheme  is; 

-  A))  -  (u"_,  +  -  A))l 

where  ^  ^  £  1  and  s"  =  m[A+u”,  A_u”].  for 

X 

m[x,y]  =  —hr(a:y)  min(lx|,  [yl) 

|ii 

where  H  is  the  Heaviside  function  . 

The  only  difference  between  (2)  and  a  conventional  second  order  accurate  ap¬ 
proximation  comes  in  the  definition  of  the  approximate  slope  s"  which  is  a  nonlinear 
function  of  ,  u" ,  . 
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The  scheme  we  are  ming  here  is  a  modified  version  of  the  scheme  developed  in 
[11].  We  use  the  pointwise  scheme  developed  by  Shu  and  Osher  [25],  The  scheme 
for  a  scalar  one  dimensional  conservation  law  is  the  following.  Let  the  conservation 
law  be  given  by: 

+  /(w)r  =  0 


We  partition  the  x-t  space  into  computational  cells  of  the  following  type: 


Then  any  conservative  scheme  can  be  written  as  [15]: 


At 


where  the  fj-^  represents  the  numerical  flxox  flowing  into  the  j-th  computational 
cell. 


The  fluxes  can  be  computed  in  the  following  fashion  as  described  in  [25],  .As¬ 
suming  that  there  is  h{Q  such  that: 

•x+^ 


1 

fiuix))  =  ^  HOdC 


then  it  follows  that: 


f{u(x))];\x  =  Xj  - 


h{x^+^)  -  h{Xj_x  ) 


Ax 


To  calculate  h(x)  let: 


-L 


Hix)  =  /  hiOdC 
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Then 


j 

KQdC,  =  Ax  ^  f{u{xj)) 

— oo 


The  values  of  H(x)  at  the  grid  points  can  be  calculated  from  the  above  and 
then  we  use  a  polynomial  interpolation  with  an  adaptive  stencil  to  reconstruct  the 
H(x)  function. We  call  the  reconstructed  polynomial  J?j^^(x).  Then  the  flux  can 
be  computed  as: 


i^j+0 


R{x)  is  constructed  in  the  following  way: 


Let  us  define  the  following  divided  differences  inductively  and  using  the  nota¬ 
tion  that  f{j)  =  f{xj) 

H'U)  =  +  1)  -  H°U))  =  fU) 

H\i)  =  + 1)  -  H'-'u))  =  -r'u) 

We  emphasize  that  there  is  no  need  to  calculate  the  value  of  H  or  any  of  its 
differences  since  everything  can  be  calculated  in  terms  of  the  divided  differences  of 
f(u).  To  choose  the  suitable  stencil  for  interpolation  we  calculate  the  Roe  speed  [21] 
according  to: 

=  /O  +  D-ZO) 

«0  +  i)-uO) 

Let  the  variable  ishift  be  the  amount  of  shift  of  the  stencil  from  the  j-th  grid  point: 


if  O.J+  ^  ^  0 

,  then  ishift  =  0, 

if  <  0 

,  then  ishift  =  1. 
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then  we  compare  the  divided  differences  and  we  move  the  stencil  according  to 


l/’(j  +  ishift)!  >  \f*{j  +  ishift  —  1)|  then  ishift=ishift  —  1 

We  continue  in  this  manner  and  let  m=j+ishift;  then  the  reconstructed  poly¬ 
nomial  can  be  written  as; 

=  H°{m)  +  H^(m){x  -  x^_^)  +  -  x^_^)(x  -  • 

The  value  of  the  flux  can  be  calculated  by  evaluating  the  first  derivative  of  R(x)  at 

X  ,  I  1  . 

The  above  scheme  is  known  to  contain  nonphysical  “expansion  shock”  solutions 
occasionally.  To  remedy  this  case  let  us  define  a  critical  cell  where  /  ( u )  has  a  zero 
in  the  interval  I  =  (uj,  U;+i].  These  are  the  only  cells  that  could  contain  expansion 
shocks. 

In  a  critical  cell  we  use  the  following  scheme: 

Let  =  max/  1/  (u)|.  Then  define  the  following  quantities  for  /  —  /  -  1  < 

s  <  j  +  r  where  r  is  the  order  of  the  scheme: 

f^i^)  =  ^(/(^)  +  (ij+Lu(s))  f~{s)  =  ^(/(.s)  -  Oj^^uis)) 

We  do  the  reconstruction  for  f~  using  ishift  =  1  and  for  using  ishiff  =  0. 

The  flux  can  be  calculated  as  a  sum  of  the  flux  moving  to  the  right  and  the 
flux  moving  to  the  left. 
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Fig  10: 

Fig  11: 

Fig  12: 

Fig  13: 

Fig  14: 

Fig  15: 

Fig  16: 

Fig  17: 

Fig  IS: 

Fig  19: 

Fig  20: 


v{t)  =  H{t),  Time  =  5,  ^  =  .2,  Ax  =  .01 
comparison  of  velocity  profile  for 
bias  voltage  1.,  1.5,  2. 
comparison  of  temperature  profile  for 
bias  voltage  1.,  1.5,  2. 

Current  vs.  Voltage 

v{t)  =  t.  Time  =  20,  ^  =  .1,  Ax  =  .01 

Evolution  of  density 

v{t)  ~  t.  Time  =  20,  ^  =  .1.  Ax  =  .01 

Evolution  of  electric  field 

v{t)  =  t,  Time  =  20,  ^  =  .1,  Ax  =  .01 

Evolution  of  temperature 

v{t)  =  t.  Time  =20.  . I,  Ax  =  .01 

Evolution  of  velocity 

v(t)  =  t.  Time  =  20,  ^  =  .1,  Ax  =  .01 

Evolution  of  Mach  number 

v{t)  =  t.  Time  =  20,  ^  =  .1,  Ax  =  .01 

Velocity  in  time 

v{t)  =  2H(t),  Time  =  3,  =  .l,Ax  =  .01.  To  =  1000/v 

evolution  of  temperature 

vit)  =  2H{t),  Time  =  3,  =  .1,  Ax  =  .01,  To  =  lOQOK 

density  in  time 
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v{t)  =  H{t),  Time  =  4,  =  .01,  Ax  =  .01,  To  =  50A' 


Fig  21;  evolution  of  mach  number 

v(t)  =  H{t),  Time  =  4,  ^  =  .01,  Ax  =  .01.  Tq  =  50A' 

Fig  22:  evolution  of  density 

v(t)  —  H{t),  Time  =  4,  ^  =  .01,  Ax  =  .01,  Tq  =  50A’ 

Fig  23:  density  in  time 

v{t)  =  2H{t),  Time  =  ^  =  .01,  Ax  =  -01, Tq  =  50A' 

Fig  24;  density 

v{t)  =  2H{t).  Time  =  1,  ^  =  .01,  Ax  =  .01,  Tq  =  50A 
Fig  25;  Current  vs.  Time 

v{t)  =  SQW{t),  Time  =  6,^  =  .1,  Ax  =  .01.  To  =  300A' 

Fig  26:  velocity  in  Time 

v{t)  =  SQW{t),  Time  =  6,  =  .1.  Ax  =  .01.  To  =  300A 

Fig  27;  Current  vs.  Time 

v(t)  =  2sm(27rt).  T ime  =  6,  ^  =  .1,  Ax  =  .01,  To  =  300A' 

Fig  2S:  density  in  time 

v{t)  —  2sm(27rt),  T ime  =  5,  ;^  =  .1,  Ax  =  .01.  To  =  300A' 

Fig  29;  temperature  in  time 

v(  t)  =  2sin(2TTt),  Time  =  b,  ^  =  .1,  Ax  =  .01.  To  =  300A' 

Fig  30:  density,  velocity,  and  momentum  not  to  scale,  third  order  stencil 
Fig  31:  density,  velocity,  and  momentum  not  to  scale,  .sixth  'xder  stencil 

Acknowledgement:  The  authors  would  like  to  thank  Dr.  Luis  Reyna  for  his 


24 


constructive  comments  on  a  first  draft  of  this  paper. 


25 


REFERENCES 


[1]  G.  Baccarani  and  M.R.  Wordeman,  “An  investigation  of  steady-state  velocity 
overshoot  effects  in  Si  and  GaAs  devices”,  Solid  State  Electr.  28  (1985),  407- 
416. 

[2]  F.J.  Blatt,  “Physics  of  electric  conduction  in  solids”,  McGraw-Hill,  New  York 
(1968). 

[3]  K.  Blotekjaer,  “Transport  equations  for  electrons  in  two-valley  semiconduc¬ 
tors”,  IEEE  Trans.  Electron  Devices  ED-17  (1970),  38-47. 

[4]  C.  Cercignani,  “The  Boltzmann  equation  and  its  application”,  Springer- Verlag, 
1988. 

[5]  D.M.  Caughey  and  R.E.  Thomas,  “Carrier  mobilities  in  silicon  empirically 
related  to  doping  and  field”,  Proc.  IEEE  55  (1967),  2192-2193. 

[6]  R.K.  Cook  and  J.  Frey,  “Two-dimensional  numerical  simulation  of  energy  trans¬ 
port  effects  in  Si  and  GaAs”,  MESFET’s,  IEEE  Trans  Electron  Devices  ED-29 
(1982),  970-977. 

[7]  R.K.  Cook  and  J.  Frey,  “An  effecient  technique  for  two-dimensional  simulation 
of  velocity  overshoot  effects  in  Si  and  GaAs  devices”,  COMPEL  1  ( 1982),  65-87. 

[8]  M.  Fischettei  and  M.  Rudan,  “Hydrodynamic  and  Monte-Carlo  simulation  of 
an  —  n  —  n"*"  device” ,  to  appear. 

[9]  C.L.  Gardner,  J.W.  Jerome,  and  D.J.  Rose,  “Numerical  methods  for  the  hydro- 
dynamic  device  model:  subsonic  flow”,  IEEE  Trans.  Computer-Aided  Design 
of  Integrated  Circuits  and  Systems,  V.  8  (1989),  501-507. 

[10]  B.  Gustafsson  and  A.  Sundstrom,  “Incompletely  parabolic  systems  in  fluid 
dynamics”,  SIAM  J.  Appl.  Math.,  V.  35,  (1978),  343-357. 

[11]  A.  Harten,  B.  Engquist,  S.  Osher  and  S.  Chakravarthy,  “Uniformly  high  order 
accurate  essentially  non-oscillatory  schemes.  III”,  J.  Comput.  Phys.,  V.  49 
(1987),  pp.231-303. 

[12]  Iv.  Hess,  “Balhstic  electron  transport  in  semiconductors",  IEEE  Trans  Electron 
Devices  ED-28  (1981),  937-  . 

[13]  K.  Hess  and  C.T.  Sah,  “Hot  electrons  in  short-gate  charge-coupled  devices". 
IEEE  Trans.  Electron  Devices,  ED-25  (1978),  1399-1405. 

[14]  S.  Laux  and  R.  Byrnes,  “Semiconductor  device  simulators  using  generalized 


26 


mobility  models”,  IBM  J.  Res.  Dev.  29  (1985),  289-301. 

[15]  P.D.  Lax  and  B.  Wendroff,  “Systems  of  conservation  laws”,  Comm.  Pure  Appl. 
Math,  V.  13,  (1960),  pp.  217-237. 

[16]  T.J.  Maloney  and  J.  Frey,  “Transient  and  steady-state  electron  transport  in 
GaAs  and  InP”,  J.  Appl.  Phys.  48  (1977),  781-787. 

[17]  C.C.  McAndrew,  K.  Singhal  and  E.L.  Heasell,  “A  consistent  nonisothermal 
extension  of  the  Scharfetter-Gummel  stable  difference  approximation”,  IEEE 
Electron  Dev.  Lett.  EDL-6  (1985),  446-447. 

[18]  J.P.  Nougier,  J.  Vaissiere,  D.  Gasquet,  J.  Zimmermarr,  and  E.  Constant,  “De¬ 
termination  of  the  transient  regime  in  semiconductor  devices  using  relaxation 
time  approximations”,  J.  Appl.  Phys.  52  (1981),  pp.  825-832. 

[19]  F.  Odeh,  M.  Rudan  and  J.  White,  “Numerical  solution  of  the  hydrodynamic 
model  for  a  one-dimensional  semiconductor  device”,  COMPEL  6  (1987),  151- 
170. 

[20]  H.D.  Rees,  “Calculation  of  distribution  functions  by  exploitng  the  stability  of 
the  steady-state”,  J.  Phys.  Chem.  Sol.  30  (1967),  pp.  643-  . 

[21]  P.  Roe,  “Approximate  Riemann  solvers,  parameter  vectors,  and  difference 
schemes”,  J.  Comput.  Phys.,  V.  43,  (1981),  pp.  357-372. 

[22]  J.  Ruch,  “Electron  dynamics  in  short  channel  field  effect  transistors”,  IEEE 
Trans.  Electron  Devices  ED-19  (1972),  pp.  652-  . 

[23]  M.  Rudan  and  F.  Odeh,  “Multi-dimensional  discretization  scheme  for  the  hy¬ 
drodynamic  model  of  semiconductor  devices”,  COMPEL  5  (1986).  149-183. 

[24]  S.  Selberherr,  “An  analysis  and  simulation  of  semiconductor  devices”,  Springer- 
Verlag,  1984. 

[25]  C.-W.  Shu  and  S.  Osher,  “Efficient  implementation  of  essentially  non-oscillatory 
shock  capturing  schemes,  II,  “J.  Comp.  Physics,  (1989),  to  appear. 

[26]  M.S.  Shur,  “Influence  of  nonuniform  field  distribution  on  frequency  limits  of 
GaAs  field  effect  transistors”.  Elect.  Lett.  12  (1976),  pp.  615-616. 

[27]  M.S.  Shur  and  E.F.  Eastman,  “Ballistic  transport  in  semiconductors  at  low 
temperatures  for  low-power  high-speed  logic”,  IEEE  Trans.  Electron  Devices 
ED-26  (1979),  pp.  1677-1683. 

[28]  R.  Stratten,  “Diffusion  of  hot  and  cold  electrons  in  semiconductor  barriers". 


27 


Phys.  Rev.  126  (1962),  pp.  2002-2014. 

[29]  T.  Tang,  “Extension  of  the  Scharfetter-Gummel  algorithm  to  the  energy  bal¬ 
ance  equation”,  IEEE  Trans.  Electron  Devices  ED-31  (1984),  pp.  1912-1914. 

[30]  K.K.  Thornber,  “Current  equations  for  velocity  overshoot”,  IEEE  Electron 
Dev.  Lett.  EDL-3  (1982),  pp.  69-71. 


28 


Density  of  Doors 


xl05 


Fig  1: 
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Fig  2:  velocity 


zero  bias,  Time  =  3,  ^  =  .2,  Ajt  =  .0 
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Time=  3 


Fig  3:  temperature 

zero  bias,  Time  =  3,  ^  =  .2,  i 
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Fig  4: 


velocity  in  steady  state 

zero  bias,  Time  =  3,  ^  =  .2,  Ax  —  .01 
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Fig  5:  n  —  riD  in  steady  state 

zero  bias,  T  =  5,  ^  =  .2,  Ax  =  .01 
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Electric  Field 


Fig  6;  electric  field  in  steady  state 

zero  bias,  Time  =  3,  ^  =  .2,  Ai  =  .01 
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ELCCrRIC  FIELD 


Timc=  3  Bias  Voltage  =  1 


Fig  8: 


electric  6eld  in  time 

„(()  =  H[i),  Time  =  5,  =  .2,  Ai  =  -01 
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VELOCITY 


Timc=  5  Bias  Voltage  =  I 


Fig  9: 


velocity  in  time 

v(t)  =  H{t),  Time  =  5,  ^  =  .2,  Ax  =  .01 
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Fig  10:  comparison  of  velocity  profile  for 


bias  voltage  1.,  l-5i  2. 
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-DENSITY 


Time=  4  Bias  Voltage  =  I 


Fig  20:  density  in  time 

t;(<)  =  H{t),  Time  =  4,  =  .01,  Az  =  .01,  To  =  50A' 
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Mach  Number 


Fig  21:  evolution  of  mach  number 

v(t)  =  Time  =  4,  =  .01,  Ai  =  .01,  To  =  50A' 
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Evolution  of  Density 
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Fig  22:  evolution  of  density 

v{i)  =  H{t),  Time  =  4,  =  .01,  Ax  =  .01,  To  =  50A' 


50 


-DENSITY 


Fig  23:  density  in  time 

v(t)  =  2H{i),  Time  =  1,  =  .01,  Ax  =  .01,  To  =  SOA" 
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micro  amp 


Current  vs.  Voltage 


Fig  12:  Current  vs.  Voltage 

„(()  =  (,  Time  =  20,  =  1,  A*  =  -01 
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Times  3  Bias  Volugc  =  2 


Fig  18:  Velocity  in  time 

v(t)  =  2H{t),  Time  =  3,  =  .1,  Ax  =  .01,  To  =  lOOOA' 
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pico  second 

Current  vs.  Time 

v{t)  =  23in(2irt),  Time 

=  6, =  .1,  Ai  =  .01,  To  =  300K 

-DENSITY 


Time=  5  Bias  Voltage  =  2 


Fig  28:  density  in  time 

t;(t)  =  2sin{2nt),  Time  =  5,^  =  .1,  Ax  =  .01,  To  =  300A“ 
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TEMPERATURE 


Time=  5  Bias  Voltage  =  2 


Fig  29:  temperature  in  time 


v(t)  =  2stn(2ff<),  Time  =  5, 


=  .l,Ax  =  .Ol.To  =300A' 


57 


Third  order  ENO 
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Fig  30:  density,  velocity,  and  momentum  not  to  scale,  third  order  stencil 
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Fig  31:  density,  velocity,  and  momentum  not  to  scale,  sixth  order  stencil 
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